Given \(n\) data points \(\{x_i, y_i\}_{i=1}^n\)

Given \(n\) data points \(\{x_i, y_i\}_{i=1}^n \xrightarrow{\text{Predict}} (x_{\texttt{new}}, y_{\texttt{new}})\)

If we have a vector of random variables \(\mathbf{x}\) and

\[ \mathbf{x} \sim \mathcal{N}_d(\boldsymbol{\mu}, \mathbf{\Sigma}) \]

then, the joint probability mass of \(\mathbf{x}\) is given by the multivariate normal:

\[ p\left( \mathbf{x} \,|\, \boldsymbol{\mu}, \mathbf{\Sigma} \right) \propto \exp \left\{ -\frac{1}{2} (\mathbf{x}-\boldsymbol{\mu})^\top \Sigma^{-1} (\mathbf{x}-\boldsymbol{\mu}) \right\} \]

A two-dimensional MVN with \(\boldsymbol{\mu}=[0,0]\) and \(\Sigma=\begin{bmatrix} 1 & 0.7\\ 0.7 & 1 \end{bmatrix}\):

# mean vector
mu <- c(0, 0)

# covariance matrix
Sigma <- matrix( c(1, 0.7, 
                   0.7, 1), nrow = 2) 

# grid of x1 and x2 values
x1x2_grid <- expand.grid(x1 = seq(-3, 3, length.out = 100), 
                         x2 = seq(-3, 3, length.out = 100))

# probability contours
probabilities <- dmvnorm(x1x2_grid, mean = mu, sigma = Sigma)

Probability contour plot:

Probability contour plot:

We can condition on one of the variables, \(p(x_2 \,|\, x_1,\,\Sigma)\)

In a simple example, imagine we are given three data points \(\{(x_i, y_i)\}_{i=1}^3\). We need to predict the function thereafter

A Gaussian process is a generalization of a multivariate Gaussian distribution to infinitely many variables.

A Gaussian process is a collection of random variables, any finite number of which have a joint Gaussian distribution.

As a distribution over functions, a Gaussian process is completely specified by two functions:

• A mean function, \(m(\mathbf{x}) = \mathbb{E}[f(\mathbf{x})]\) and
• A covariance function, \(k(\mathbf{x}, \mathbf{x^\star}) = \mathbb{E}\Big[\big(f(\mathbf{x})- m(\mathbf{x})\big) - \big(f(\mathbf{x^\star})- m(\mathbf{x^\star}\big) \Big]\)

\[ f(\mathbf{x}) \sim \mathcal{GP}\big(\, m(\mathbf{x}), k(\mathbf{x}, \mathbf{x^\star}) \, \big) \]

Generative model

\[ y(\mathbf{x}) = f(\mathbf{x}) \Big[ + \epsilon\sigma_y \Big]\\ p(\epsilon)=\mathcal{N}(0,1) \]

Place GP prior over the nonlinear function (mean function often taken as 0).

\[ p(f(\mathbf{x}) \, | \, \theta) = \mathcal{GP}\big(0, k(\mathbf{x}, \mathbf{x^\star})\big)\\ k(\mathbf{x}, \mathbf{x^\star}) = \sigma^2 \exp \left\{ -\frac{1}{2\ell^2}(x-x^\star)^2 \right\} \]

\[ p(y_1, y_2) = \mathcal{N} \left( \begin{bmatrix} \boldsymbol{\mu_1\\ \mu_2} \end{bmatrix}, \begin{bmatrix} \mathbf{\Sigma}_{11} & \mathbf{\Sigma}_{12}\\ \mathbf{\Sigma}_{21} & \mathbf{\Sigma}_{22} \end{bmatrix}\right) \] Bayes’ theorem:

\[ p(\mathbf{y}_1 \,|\, \mathbf{y}_2) = \frac{p(\mathbf{y}_1, \mathbf{y}_2)}{p(\mathbf{y}_2)} \]

With an involved proof:

• Predictive mean (linear projection in inverse distance space): \(\boldsymbol{\mu}_{\mathbf{y}_1|\mathbf{y}_2} = \boldsymbol{\mu_1} + \mathbf{\Sigma}_{12}\mathbf{\Sigma}_{22}^{-1}(\mathbf{y}_2-\boldsymbol{\mu}_2)\)
• Predictive variance (quadratic reduction in variance): \(\mathbf{\Sigma}_{\mathbf{y}_1|\mathbf{y}_2} = \mathbf{\Sigma}_{11} - \mathbf{\Sigma}_{12}\mathbf{\Sigma}_{22}^{-1}\mathbf{\Sigma}_{21}\)

• Premise: We have an unknown and expensive-to-evaluate objective function, \(y(x)\)
• Goal: Return an value \(x_M\) which corresponds to a minimum \(y(x_M)\)
• High evaluation cost requires selecting the location of new evaluations carefully

In Bayesian optimization, a surrogate function (usually called an acquisition function) is used to select the next best evaluation location

Steps:

1. Evaluate the objective function at a point.
   
2. Fit/update Gaussian process model with the data.
   
3. Optimize an acquisition function to select the next sampling point.
   
4. Repeat until a defined convergence or exhaustion of resources.

• We are minimizing our objective function \(y(x)\).
• Call the best solution so far \(x_M^\star\).
• Improvement can be defined as:

\[ I(x_{\texttt{next}}) = \max\left\{0, f(x_M^\star) - \widetilde{f}(x_{\texttt{next}}) \right\} \]

Each point \(x_{\texttt{next}}\) in a Gaussian process is associated with a mean given by evaluating the mean function at \(x_{\texttt{next}}\), \(\mu(x_{\texttt{next}})\) and a variance given by evaluating the variance function at \(x_{\texttt{next}}\), \(k(x, x_{\texttt{next}})\).

• Function evaluations so far: \((\mathbf{x}_0, \mathbf{y}_0)\)
• Define \(\eta := \min \mathbf{y}_0\)
• Define loss function:

\[ \lambda(y) := \begin{cases} y;\quad y < \eta \\ \eta;\quad y \ge \eta \end{cases} \]

Fancy way of saying loss is \(\min(y, \eta)\).

\[ \mathbb{E}[\lambda(y)] := \int_{y=-\infty}^{y=\infty} \lambda(y)\cdot p(y \,| \,x, I_0) \, dy \]

where \(x\) is the new selected point and \(I_0\) is the current “knowledge” (e.g., Gaussian process parameters).

Our integral must be split at the critical value \(\eta\). When \(y<\eta\), \(\lambda(y)=y\) and when \(y\ge\eta\), \(\lambda(y)=\eta\).

\[ \mathbb{E}[\lambda(y)] := \int_{y=-\infty}^{y=\eta} y\cdot p(y \,| \,x, I_0) \, dy + \int_{y=\eta}^{y=\infty} \eta\cdot p(y \,| \,x, I_0) \, dy \]

Note that \(p(y \,| \,x, I_0)\) is driven by the Gaussian process!

\[ p(y \,| \,x, I_0) = \mathcal{N}\big(y; m(y\,|\,I_0), k(y_0,y\,|\,I_0)\big) \]

Plugging this into the integrals, we get

\[ \begin{multline} \mathbb{E}[\lambda(y)] = \int_{-\infty}^{\eta} y\cdot \mathcal{N}\big(y; m(y\,|\,I_0), k(y_0,y\,|\,I_0)\big) \, dy \\ + \int_{\eta}^{\infty} \eta\cdot \mathcal{N}\big(y; m(y\,|\,I_0), k(y_0,y\,|\,I_0)\big) \, dy \end{multline} \]

\[ \begin{multline} =\eta + (m(y\,|\,I_0) - \eta) \Phi(\eta\,|\,m(y\,|\,I_0),k(y_0,y\,|\,I_0)) \\ - k(y_0,y\,|\,I_0)\phi(\eta\,|\,m(y\,|\,I_0),k(y_0,y\,|\,I_0) ) \end{multline} \]

• Strongly correlated observations makes for poorly conditioned covariance matrices
• Derivative observations can be seamlessly incorporated into a GP
• A function over which we have a GP is jointly Gaussian with its derivatives
• Derivative observations can be treated just as observations of the function itself (with modified covariance function)
• A derivative is only weakly correlated with the function very nearby, improving conditioning

GPGO was used to find the optima for several known, difficult functions.

Performances was measured using “gap”:

\[ G := \frac{y(x^{\texttt{first}}) - y(x^{\texttt{best}})}{y(x^{\texttt{first}}) - y^{\texttt{opt}}} \]

• \(y^{\texttt{opt}}\) is the known global optimum
• \(x^{\texttt{first}}\) and \(x^{\texttt{best}}\) are the first point evaluated and best point found

• Work does not explore optimization problems with which the new GPGO method might struggle or fail
• Explore reasons for lower performance of GPGO compared to other methods when optimizing certain functions such as Shu or A5
• Plan to apply BO methods to computationally intensive objection functions in clinical trial design

Journal article: Osborne et al. (2009)

Mathematics/derivations: Williams and Rasmussen (2006), Osborne et al. (2009)

Inspiration for some of the visualizations: Turner (2016), Gramacy (2022a), and Gramacy (2022b)

Bayesian optimization visualization: Farina (2024)